Shear modulus of simulated glass-forming model systems: 
Effects of boundary condition, temperature and sampling time 
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The shear modulus G of two glass-forming colloidal model systems in d = 3 and d = 2 dimensions is 
investigated by means of, respectively, molecular dynamics and Monte Carlo simulations. Comparing 
ensembles where either the shear strain 7 or the conjugated (mean) shear stress r are imposed, we 
compute G from the respective stress and strain fluctuations as a function of temperature T while 
keeping a constant normal pressure P. The choice of the ensemble is seen to be highly relevant for the 
shear stress fluctuations I-If{T) which at constant r decay monotonously with T following the affine 
shear elasticity ij,a{T), i.e. a simple two-point correlation function. At variance, non-monotonous 
behavior with a maximum at the glass transition temperature Tg is demonstrated for /iF(r) at 
constant 7. The increase of G below Tg is reasonably fitted for both models by a continuous cusp 
singularity, G{T) cx (1 — T/Tg)^^'^ , in qualitative agreement with some recent replica calculations. 
It is argued, however, that longer sampling times may lead to a sharper transition. The additive 
jump discontinuity predicted by mode-coupling theory and other replica calculations thus cannot 
ultimately be ruled out. 

PACS numbers: 61.20. Ja,65.20.-w 



I. INTRODUCTION 

Stress fluctuation formalism. Among the fundamen- 
tal properties of any solid material are its elastic con- 
stants [ll-Q- They describe the macroscopic response 
to weak external deformations and encode information 
about the potential landscape of the system. In their 
seminal work Q Squire, Holt and Hoover derived an ex- 
pression for the isothermal elastic constants extending 
the classical Born theory [l| to canonical ensembles of 
imposed particle number TV, volume V and finite (mean) 
temperature T. They found a correction term Cp^^^ to 
the Born expression for the elastic constants which in- 
volves the mean-square fluctuations of the stress 0-01 
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with l3 — l/k^T being the inverse temperature, /cb Boltz- 
mann's constant, aap the instantaneous value of a com- 
ponent of the stress tensor cTap = (o'ap) and Greek letters 
denoting the spatial coordinates in d dimension. Perhaps 
due to the fact that the demonstration of the "stress 
fluctuation formalism" in Ref. Q assumed explicitly a 
well-defined reference position and a displacement field 
for the particles, it has not been appreciated that this 
approach is consistent (at least if the initial pressure is 
properly taken into account) with the well-known pres- 
sure fluctuation formula for the compression modulus K 
of isotropic fluids [6-9] which has been derived several 
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decades earlier, presumably for the first time in the late 
1940s by Rowlinson In any case, the stress fiuctuation 
formalism provides a convenient route to calculating the 
elastic properties in cornputer simulations which has been 
widely used in the past j?!. [l0l - l2l| . It has been generalized 
to systems with nonzero initial stress [l3| , hard-sphere in- 
teractions [l3] and arbitrary continuous potentials [sj , or 
to the calculation of local mechanical properties [20] • 

Low-temperature limit and non-affine displacements. 
In particular from Ref. 5] it has become clear that the 
stress fiuctuation correction to the leading Born term 
does not necessarily vanish in the zero-temperature limit. 
This is due to the fact that when a system is sub- 
jected to a homogeneous deformation, the ensuing parti- 
cle displacements need not follow the macroscopic strain 
affinely [l3l - [T6l [T8| . The stress fluctuation term quanti- 
fies the extent of these non-affine displacements, whereas 
the Born term refiects the affine part of the particle dis- 
placements. How important the non-affine motions are, 
depends on the system under consideration [Tsj. While 
the elastic properties of crystals with one atom per unit 
cell are given by the Born term only, stress fluctuations 
are significant for crystals with more complex unit cells 
[^]. They become particularly pronounced for poly mer 
like soft materials [2l| and amorphous solids fH. Illl [isl - 
[Tsl . [20I [22} . This is revealed by recent simulation stud- 
ies of various glass formers [ll| , like Lennard- Jones (L J) 
mixtures jisl-llq] , polymer glasses 0, HOl or silica melts 
i22j. Many of these simulation studies concentrate on 
the mechanical properties deep in the glassy state, ex- 
ploring for instance correlations between the non-affine 
displacement field and vibrational anomalies of the glass 
( "Boson peak" ) [l^ , the mechanical heterogeneity at 
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the nanoscale [20|, |23| or the onset of moleeular plasticity 
in the regime where the macroscopic deformation is still 
elastic [24| . 

Shear stress and shear modulus. Surprisingly, the 
temperature dependence of the elastic constants of glassy 
materials does not appear to have been investigated ex- 
tensively using the stress fluctuation formalism [3, ll.[l9|. 
Following the pioneering work by Barrat et al. we 
present in this work such a study where we focus on the 
behavior of the shear modulus G of two isotropic colloidal 
glass-formers in d = 3 and d = 2 dimensions sampled, re- 
spectively, by means of molecular dynamics (MD) and 
Monte Carlo (MC) simulations 0,0, [HI. Since we are 
interested in isotropic systems one can avoid the incon- 
venient tensorial notation of the full stress fluctuation 
formalism as shown in Fig. [TJ Focusing on the shear in 
the xy-plane we use r = a^y for the mean shear stress 
and T = dxy for its instantaneous value. As an important 
contribution to the shear modulus we shall monitor the 
shear stress fluctuations 



(2) 



as a function of temperature T. Quite generally, the shear 
modulus G may then be computed according to the stress 
fluctuation formula 0, 



Gtt = MA — Mf 



(3) 



where /ia stands for the so-called "afRne shear elasticity" 
corresponding to an assumed ajjine response to an exter- 
nal shear strain 7. For pairwise interaction potentials it 
can be shown that 



MA = Mb - -Pcx > 0, 



(4) 



i.e. MA comprises the so-called "Born-Lame coefhcient" 
MB) a simple moment of the first and the second deriva- 
tive of the pair potential with respect to the distance of 
the interacting particles [11, 0, 13, and the excess con- 
tribution Pox to the total pressure P. These notations 
will be properly defined in Sec. Ill B I below where a sim- 
ple derivation will be given. Since Eq. ^ is not the only 
operational definition of the shear modulus discussed in 
this work, an index has been added reminding that the 
shear stress f is used for the determination of G. 

Physical motivation and systems of interest. If an 
infinitesimal shear stress r is applied at time t — 0, 
an isotropic solid exhibits, quite generally, an elastic 
response quantified in strength by a finite static {t- 
independent) shear modulus (G > 0) as shown by the 
top curve in Fig. [IJb), whereas a liquid has a vanish- 
ing modulus (G — 0) and flows on long time scales as 
shown by the bottom curve. The shear modulus G can 
thus serve as an "order parameter" distinguishing liq- 
uid and solid states [26i - i31i] . Upon melting, the shear 
modulus of crystalline solids is known to vanish discon- 
tinuously with T. A natural question which arises is that 
of the behavior of G{T) for amorphous solids and glasses 
in the vicinity of the glass transition temperature Tg. 
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FIG. 1: Setup, notations and ensembles investigated numeri- 
cally: (a) Schematic experimental setup for probing the shear 
modulus G. (b) Expected shear modulus G{t) as a function of 
measurement time t ai T <^Tg (top curve), T slightly below 
Tg (middle curve) and in the liquid limit for T ^ (bottom 
curve), (c) We use periodic simulation boxes with deformable 
box shapes where we impose r = 0. The shear modulus G 
may be computed using Eq. © or Eq. from the instanta- 
neous strains 7 and stresses f . (d) Using Eq. Q G may also 
be computed using the afSne shear elasticity /iA and the shear 
stress fluctuations /ip in ensembles of fixed strain (7 = 0). (e) 
Combined volumetric and shear strain fiuctuations allowing 
the computation of both compression modulus K and shear 
modulus G from the respective strain fluctuations. 



Interestingly, qualitative different theoretical predictions 
have been put forward by mode-coupling th eory (MCT) 
[31I . [3^ and some versions of replica theory [sO] predict- 
ing a discontinuous jump at the glas s transition whereas 
other versions of replica theory |28l . [29| suggest a con- 
tinuous transition in the static limit of large sampling 
times t. (See Ref. [29[ for a topical recent discussion of 
the replica approach.) The numerical characterization of 
G(r) for colloidal glass-former is thus of high current in- 
terest. Since the stress fluctuation formalism can readily 
be added to the standard simulation codes of colloidal 
glass-formers [ll|j [Hj 113 imposed particle number TV, 
box volume V, box shape (7 — 0) and temperature T, 
called NV7T-ensembles, this suggests the use of Eq. ([3]) 
to tackle this issue. Note that the latter stress fluctuation 
relation should also be useful for wide range of related 
systems, beyond the scope of the present work, including 
the glass transition of polymeric liquids [7] , colloidal gels 
[S^ l and self-assembled networks [5^ as observed exper- 
imentally fo r hy perbranched polymer chains with sticky 
end-groups [35] or bridged networks of telechelic poly- 
mers in water-oil emulsions 1361 1371 . 
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Validity through a solid-liquid transition. This begs 
the dehcate question of whether the general stress fluc- 
tuation formalism and especially Eq. ([3]) for the shear 
modulus only holds for elastic solids with a well-defined 
displacement field, as suggested in the early literature j3- 
0] , or if the approach may actually be also applied through 
a solid-liquid transition up to the liquid state [1, [2l| with 
MA and thus G — Grr —J' as it should for a liq- 
uid. By reworking the theoretical arguments leading to 
Eq. Q it will be seen that this is in principal the case, 
at least if standard thermodynamics can be assumed to 
apply. The latter point is indeed not self-evident for the 
colloidal glasses we are interested in where (i) some de- 
grees of freedom are frozen on the time scale probed and 
(ii) the measured response to an imposed infinitesimal 
shear t may depend on the sampling time t, especially 
for temperatures T around Tg as shown in Fig. [Ijb). To 
be on the safe side, it thus appears to be appropriate 
to test numerically the applicability of the stress fluc- 
tuation formula, Eq. ([3]), using the corresponding strain 
fluctuation relations which are more fundamental on ex- 
perimental grounds. 

Strain fluctuation relations. Panel (a) of Fig. [1] shows 
a schematic setup for probing experimentally the shear 
modulus G of isotropic systems confined between two 
rigid walls. We assume that the bottom wall is fixed and 
the particle number N, the (mean) normal pressure P 
and the (mean) temperature T are kept constant. In lin- 
ear response either a shear strain 7 — tan(Q;) or a (mean) 
shear stress r may be imposed. The thin solid line indi- 
cates the original unstrained reference system with 7 = 
and r = 0, the dashed line the deformed state. From both 
the mechanical and the thermodynamical point of view, 
the shear modulus is defined as 



G 



(5) 



7^0 



where the derivative may also be taken at r = 0. In 
practice, one obtains G by fitting a linear slope to the 
(7, r)-data at vanishing strain. Obviously, the noise level 
is normally large for small strains. From the numerical 
point of view it would thus be more appropriate to work 
at imposed mean stress t = 0, to let the strain freely 
fluctuate and to sample the instantaneous values (7, t). 
The shear modulus G is then simply obtained by linear 
regression 



G. 



(SjSt) 



(572 



(6) 



T = 



where the index indicates that both the instantaneous 
values 7 and f are used. (Note that imposing the strain 
7 also fixes the instantaneous strain 7 which thus can- 
not fluctuate.) The associated dimensionless correlation 
coefficient 



(7) 



allows to determine the quality of the fit. We do not 
know whether Eq. ^ and Eq. (O have actually been 
used in a real experiment, however since no difficulty 
arises to work in a computer simulation at constant r 
and to record (7,t), we use this linear regression as our 
key operational definition of G which is used to judge 
the correctness of the stress fluctuation relation, Eq. ([3]) , 
computed at constant 7. To avoid additional physics at 
the walls, we use of course periodic boundary conditions 
9] (the central image being sketched) as shown by the 
panels on the right-hand side of Fig. [TJ As discussed in 
Sec. Ill C| it is possible to determine G^r at constant vol- 
ume (NVrT-ensemble) as shown by panel (c) or constant 
normal pressure (NPrT-ensemble) as shown by panel (e). 
The latter ensemble has the advantage that one can eas- 
ily impose in addition the same normal pressure P for all 
temperatures T as it is justified for experimental reasons. 
If thermodynamics holds, it is known (Sec. Ill Ap that 



{5^5f)\^^kBT/V 



(8) 



should hold [Sgl- Using Eq. ^ this implies a second 
operational definition 



G 



_ kBT/V 



77 



(57^) 



(9) 



which has the advantage that only the instantaneous 
shear strain 7 has to be recorded. The main technical 
point put forward in this paper is that all operational 
definitions yield similar results 



G 



7-rlr=0 



G 



77It=0 



7=0 ' 



(10) 



even for temperatures where a strong dependence on the 
sampling time t is observed. In other words, Eq. (|10p 
states that this time dependence applies similarly (albeit 
perhaps not exactly) to the various moments computed 
for the three operational definitions of G in different en- 
sembles. 

Outline. The paper is organized as follows: Several 
theoretical issues are regrouped in Sec. HT] We begin the 
discussion in Sec. Ill Al by reminding a few useful gen- 
eral thermodynamic relations. A simple derivation of the 
stress fluctuation formula, Eq. ([3|), for the shear modu- 
lus Grr in constant-7 ensembles is given in Sec. IIIBI 
(The analogous stress fluctuation formula for the com- 
pression modulus K in constant- ensembles @ is red- 
erived in Appendix lA 21 ) Shear stress fluctuations /xp 
in ensembles with different boundary conditions are dis- 
cussed in Sec. Ill Cl For constant-r ensembles /ip = /iFl,- 
will be shown to reduce to the affine shear elasticity /.tA, 
i.e. the functional Grr must vanish for all temperatures. 
As pointed out recently 39] , the truncation of the inter- 
action potentials used in computational studies implies 
an impulsive correction for the affine shear elasticity /xa. 
This is summarized in Sec. Ill Dl The numerical models 
used in this study are presented in Sec. lIIII We turn then 
in Sec. IIVI to our main numerical results. Starting with 
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the high-T liquid limit, Sec. lIVAl shows that all our oper- 
ational definitions correctly yield a vanishing shear mod- 
ulus. Moving to the opposite low-T limit we present in 
Sec. IIVBI the elastic response for a topologically fixed 
network of harmonic springs constructed from the "dy- 
namical matrix" of our 2D model glass system [l^ [ij] . 
As shown in Sec. IIV CI the same behavior is found quali- 
tatively for low-T glassy bead systems. The temperature 
dependence of stress fluctuations and shear moduli in dif- 
ferent ensembles is then discussed in Sec. IIV D] As stated 
above, one central goal of this work is to determine the 
behavior of the shear modulus around Tg. We conclude 
the paper in Sec. [V] Thermodynamic relations and nu- 
merical findings related to the compression modulus K 
are regrouped in Appendix |Al 



with j3 = being the inverse temperature. For our 

variables X — and I = t, Eq. p5|) implies Eq. (|8]) 
from which it is seen that Eq. ^ is consistent with 
Eq. © . The latter relation is also obtained directly using 
Eq. (HH). 

As discussed in textbooks HO] a simple average 
A = {A) of some observable A does not depend of 
the chosen ensemble, at least not if the system is large 
enough {N,V — >■ oo). A correlation function {SA6B) of 
two observables A and B may differ, however, depending 
on whether X or I are imposed. As demonstrated by 
Lebowitz, Percus and Verlet in 1967 one verifies that 



SASBj 



d{(31) dA dB 
dX d{(3I)d{pi)' 



(17) 



II. THEORETICAL CONSIDERATIONS 

A. Reminder of thermodynamic relations 

We assume in the following that standard thermody- 
namics and thcrmostatistics [40| can be applied and re- 
mind a few useful relations. Let us denote an exten- 
sive variable by X and its instantaneous value by X, the 
conjugated intensive variable by / and its instantaneous 
value by J [4l|. With F(T,X) being the free energy at 
temperature T and imposed A, the intensive variable / 
and the associated modulus M are given by 



9F(T, X) 
dX ' 



dl d^F(T,X) 

M = V = V — — - 

~ dX dX^ 



(11) 

(12) 



Alternatively, one may consider the Legendre transfor- 
mation iJ(T, /) = F(T, X) - XI for which 



X 

V 

M 



dH{T, I) 



dl 



dX 
'dl 



(13) 
(14) 



In our case, the extensive variable is A = V"f, the inten- 
sive variable I = t [42] and the modulus M — G [43| . 

We emphasize that there is an important difference be- 
tween extensive and intensive variables: If the extensive 
variable X is fixed, X cannot fluctuate. If, on the other 
hand, the intensive variable / = (/) is imposed, not only 
X but also I may fluctuate [4^]. In our case this means, 
that 7 = 7 for the NV7T-ensemble, while the shear stress 
fluctuations /ip defined above do clearly not vanish in the 
NVrT-ensemble shown in Fig. [Ud). Assuming / to be 
imposed, it is well known that [40| 



6X{l3SI)j 



SX' 



1 (15) 
(16) 



9(/3/) (3M 



If the extensive variable and at least one of the observ- 
ables are identical, the left-hand side of Eq. p7|) vanishes: 
If A = X and B = /, this yields Eq. (US]). UA = B = X, 
this implies Eq. (ITBl) . More importantly, for A = B = I 
one obtains using Eq. (|12l) 



M = PV {SI 



-I3V (SI 



X 



(18) 



If it is possible to probe the intensive variable fluctua- 
tions in both ensembles, this does in principle allow to 
determine the modulus M. (The latter equation might 
actually be taken as the definition of M in cases where 
the thermodynamic reasoning becomes dubious.) We em- 
phasize that the modulus M is an intrinsic property of 
the system and does not depend on the ensemble, albeit it 
is determined using ensemble depending correlation func- 
tions. For a thermodynamically stable system, M > 0, 
Eq. HH]) implies that 



SF 



> Si 



(19) 



and that both fluctuations must become similar in the 
limit where the modulus M becomes small. For the shear 
stress fluctuations iip, Eq. ([T8|) corresponds to the trans- 
formation 



(20) 



and one thus expects /ip to be larger for NVrT-systems 
than for NV7T-systems, while in the liquid limit where 
G — > the imposed boundary conditions should not mat- 
ter. We shall verify Eq. (|20l) numerically below. A glance 
at the stress fluctuation formula, Eq. (|3]), suggests the 
identity 



MfI^ = MA- 



(21) 



Assuming Eq. (jH]) and Eq. (|2T]) . it follows for the corre- 
lation coefficient c-yr defined above that 



(22) 
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If the measured (7, f ) were indeed perfectly correlated, 
i.e. c-yr = 1, this implies G = fj,A and thus ^J-F\v■y — 0- 



In fact, for all situations studied here we always have 



< 1, thus G < /iA, 



(23) 



i.e. the afhne shear elasticity sets only an upper bound to 
the shear modulus and a theory which only contains the 
afhne strain response must overpredict G. We show now 
directly that Eq. © and, hence, Eqs. (I21I22I23|) indeed 
hold and this under quite general conditions. 



B. Shear modulus at imposed shear strain 

Introduction. As shown in Fig.[IJa), we demonstrate 
Eq. ([3]) from the free energy change associated with an 
imposed arbitrarily small pure shear strain 7 in the xy- 
plane assuming a constant particle number N , a con- 
stant volume V and a constant temperature T (NV7T- 
ensemble). It is supposed that the total system Hamilto- 
nian may be written as the sum of a kinetic energy and 
a potential energy. Since for a plain shear strain at con- 
stant volume the ideal free energy contribution does not 
change, i.e. is irrelevant for r and G, we may focus on 
the excess free energy contribution 



Fex(r,7) = -fcBTln(Zex(7)) 



(24) 



due to the conservative interaction energy of the parti- 
cles. The (excess) partition function Zox(O) of the un- 
perturbed system at 7 = is the Boltzmann-weighted 
sum over all states s of the system which are accessible 
within the measurement time t. The argument (0) is a 
short-hand for the unstrained reference. Following the 
derivation of the compression modulus K by Rowlinson 
Q the partition function 



^cx(7) = Ee^P(-/^t^'^(^)) 



(25) 



of the sheared system is supposed to be the sum over the 
same states s, but with a different metric [46| correspond- 
ing to the macroscopic strain which changes the total in- 
teraction energy Us{'y) of state s and, hence, the weight 
of the sheared configuration for the averages computed. 
This is the central hypothesis made. Interestingly, it is 
not necessary to specify explicitly the states of the unper- 
turbed or perturbed system, e.g., it is irrelevant whether 
the particles are distinguishable or not or whether they 
have a well-defined reference position for defining a dis- 
placement field. 

Shear stress and modulus for general potential. As- 
suming Eq. (j24[) to hold we compute now the mean shear 
stress T and the shear modulus G for a general interaction 
potential Usij). We note for later convenience that 



51n(Ze.(7)) ^L(7) 



a^ln(Zex(7)) 



Zcx(7) V^cx(7) 



(26) 
(27) 



for the derivatives of the free energy and 

9^ex(7) _ ^ aTT'/..\ ^-PUAl) 



972 



-Y^PU'A^) e 

S 

S 



(28) 



(29) 



for the derivatives of the excess partition function where 
a prime denotes the derivative of a function f{x) with 
respect to its argument x. Using Eq. and taking 
finally the limit 7 — >■ one verifies for the shear stress 
that 



1 



= (f) withf^- [7^(7) 1^.0 



(30) 



defining the instantaneous shear stress Note that 

the average taken is defined as 



1 



^cx(O) 



(31) 



using the weights of the unperturbed system. The shear 
stress thus measures the average change of the total inter- 
action energy Us{^) taken at 7 = 0. The shear modulus 
G is obtained using in addition Eq. (f27| and Eq. ([29|) 
and taking finally the 7—^0 limit. Confirming thus the 
stress fluctuation formula Eq. ([3]) stated in the Introduc- 
tion, this yields 



G — Gtt = ^J-A — f-F 



7=0 



with /iA being the "afhne shear elasticity" 



MA 



1 

V 



7=0 



(32) 



(33) 



already mentioned above and /ip as defined in Eq. 
The comparison of Eq. (15^ and Eq. (PH)) confirms 
Eq. (EI]). 

Comment. The affine shear elasticity /ia corresponds 
to the change (second derivative) of the total energy 
which would be obtained if one actually strains affinely 
in a computer simulation a given state s without allow- 
ing the particles to relax their position. As shown for 
athermal (T 0) amorphous bodies [H, [ij, [l^ , the po- 
sitions of the particles of such a strained configuration 
will of course in general change slightly to minimize the 
interaction energy relaxing thus the elastic moduli. This 
is also of relevance for thermalized solids where the non- 
affine displacements of the particles are driven by the 
minimization of the free energy. It is for this reason that 
the shear-stress fluctuation term > must occur in 
Eq. (I32p correcting the overprediction of the shear mod- 
ulus G by the affine shear elasticity ha, Eq. (l23l) . This 
point has been overlooked in the early literature Jj] and 
only appreciated much later d, d, H, [l^, [3 lU as dis- 
cussed in Barrat's review [l^. Interestingly, as has been 
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shown by Lutsko Q , and other similarly defined stress 
fluctuations become temperature independent in the har- 
monic ground state approximation for T — > 0. Probing 
the stress fluctuations in a low-temperature simulation 
allows thus to determine the elastic moduli of athermal 
sohds. 

Please note that up to now we have not used explicitly 
the coordinate transformations (metric change) associ- 
ated with an affine shear strain. This is needed to obtain 
operational definitions for the instantaneous shear stress 
f (and thus for t and /ip) and the affine shear elasticity. 

Coordinate transformation. As shown in Fig. [ija), 
we assume a pure shear strain which transforms the x- 
coordinate of a particle position or the distance between 
two particles as 



x{0) x{-f) = x{0)+-fy{0) 



(34) 



leaving all other coordinates unchanged [46]. The 
squared distance between two particles thus trans- 
forms as 

r^{0) r^{"f) = r^{0) + 2jx{0)y{0) + j^y^{0) (35) 

where we need to keep the 7^-term for calculating the 
shear modulus below. We note for later reference that 
this implies 



= 2a;(0)2/(0) + 27y2(0) (36) 
= 2y2(0) (37) 



for, respectively, the first and the second derivative with 
respect to 7. Let us now consider an arbitrary function 
/(r(7)). Using Eq. (|36)) it is readily seen that its first 
derivative with respect to 7 may be written as 

^^^IP = /'(r)^(22;(0)y(0) + 27y2(0)) (38) 



rf'ir) n-j-riy for 7 — > 0. 



(39) 



In the second step we have introduced the components 
Tlx = a;(0)/r(0) and riy — j/(0)/r(0) of the normalized 
distance vector between both particles. More generally, 
we denote by ria the a-component of a normalized vector 
with a — x,y^ . . . Stating only the small-7 limit for the 
second derivative of f{r{'^)) with respect to 7 we note 
finally that 



7->0 \ C*7'' 



) = {r'nr)~rnr)) 

+ r/'(r) UyUy (40) 



2 2 

n^riy 



where we have dropped the argument (0) for the distance 
r of the unperturbed reference system. 



Pair interaction potentials. We assume now that the 
interaction energy Us^j) of a configuration s is given by a 
pair interaction potential u{r) between the particles [al^ 
with r(7) being the distance between two particles, i.e. 



Ush) = J2u{riij)) 
I 



(41) 



where the index I labels the interaction between the par- 
ticles i and j with i < j. Using the general result Eq. ([50)1 
for the shear stress it follows using Eq. ([55)) that 



n) nx,i riy^i 



(42) 



where we have dropped the argument (0) in the 7^0 
limit. We have thus rederived for the shear stress (a — 
X, (3 — y) the well-known Kirkwood expression for the 
general excess stress tensor Q 

<p = i^Tfi) with a^"^ riu'i'^ihc,,in^,i- (43) 



The affine shear elasticity /ia is obtained using Eq. ((33 
and Eq. (gn]). This yields 



MA = Mb + cr. 



vv 



(44) 



with the first contribution /iB being the so-called "Born- 
Lame coefficient" \M 



Mb 




rfu"{ri) - riu'iri)) nljul i 



(45) 



It corresponds to the first term in Eq. (PH)) . The second 
contribution cr™ stands for the normal excess stress in 
the y-direction, i.e. the a = j3 = y component of the 
excess stress tensor. 

Isotropic systems. In this work we focus on isotropic 
materials under isotropic external loads. The normal (ex- 
cess) stresses aj^^ for all directions a are thus identical, 
i.e. 



^nu'{ri) 



(46) 



with Pcx = P — Pid being the excess pressure, P the 
total pressure, Pid = k^Tp the ideal gas pressure and d 
the spatial dimension. This implies that the affine shear 
elasticity is given by ma = Mb — -Pex as already stated 
in the Introduction, Eq. We note that for isotropic 
d-dimensional systems it can be shown [21I [39| that the 
Born-Lame coefficient can be simplified as 



1 



Mb 



d{d + 2)V 



yir^u"{r,) -nu'{r,) 



(47) 



The d-dependent prefactor stems from the assumed 
isotropy of the system and the mathematical formula [i^ 



(48) 
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{6ap being the Kronecker symbol [48]) for the compo- 
nents of a unit vector in d dimensions pointing into arbi- 
trary directions. By comparing the excess pressure Pcxj 
Eq. (|46l) , with the underhned second contribution to the 
Born coefficient in Eq. (|Tf)) this imphes that 



-1 



d{d + 2)V 



Y^nu'in)) =P,^/{d + 2). 



(49) 



It is thus inconsistent to neglect the explicit excess pres- 
sure in Eq. Q while keeping the second term of the Born- 
Lame coefficient This approximation is only justified 
if the excess pressure Pcx is negligible (and not the total 
pressure P). 



C. Shear stress fluctuations in diflFerent ensembles 

Different ensembles. Being simple averages neither 
Pcx, /^A or /iB does depend for sufficiently large systems 
on the chosen ensemble, i.e. irrespective of whether 7 or 
T is imposed one expects to obtain the same values. As 
already reminded in Sec. Ill Al this is different for fiuctua- 
tions in general [9]. For instance, it is obviously pointless 
to use the shear-strain fluctuation formulae G^t or G^^ 
in a constant-7 ensemble. A more interesting result is 
predicted if Grr = MA — /^f is computed in the "wrong" 
constant-T ensemble. Since /xf|^ = Ma, Eq. (^01) . one ac- 
tually expects to observe Grr — for all temperatures T 
if our thermodynamic reasoning holds through the glass 
transition up to the liquid state. We will test numerically 
this non-trivial prediction in Sec. IIVI 

Up to now we have only considered for simplicity of 
the presentation NV7T- and NVrT-ensembles at con- 
stant volume V . Since on general grounds pure deviatoric 
and dilatational (volumetric) strains are decoupled (both 
strains commute) , one expects the observable Grr also to 
be applicable in the NP7T-ensemble shown in Fig. [Ud) 
and the observables G^r and G^^ to be applicable in the 
NPrT-ensemble as sketched in Fig.[TJe). In the latter en- 
semble one also expects to find Grr = and ii-p\r = Ma- 
This will also be tested below. 

Direct derivation of /ip for imposed t . We demon- 
strate now directly Eq. ([2T|) for the shear stress fiuctua- 
tions using the general mathematical identity 



^f{x)A[x))=-k^T{^^[x)) 



(50) 



with X being an unconstrained coordinate, A(x) some 
property, fix) = —u'{x) a "force" with respect to 
some "energy" u{x) and the average being Boltzmann 
weighted, i.e. (...) oc / da; . . . e~^'^'^^\ (It is assumed in 
Eq. ([50]) that u{x) decays sufficiently fast at the bound- 
aries.) Following work by Zwanzig j49| we express first 
the instantaneous shear stress f given above by the Kirk- 
wood formula, Eq. (I42p . by the alternative virial repre- 
sentation 

1 



Vifx,: 



(51) 



The second term, sometimes called the "inner virial", 
stands for a sum over all particles i with yi being the 
y-coordinate of the particle and f^^i the x-coordinate of 
the force acting on the particle. This contribution indeed 
vanishes on average as can be seen using Eq. (j50p . 



{yifx-: 



0, 



(52) 



which shows that (f ) = r as it must. The stress fluctua- 
tions may thus be expressed as 



MF = PV{{T-Tf) 



dyjfx.„ 
dxi 



(53) 



where we have used the integration by part, Eq. (j51l) . 
This step requires that all particle positions are un- 
constrained and independent (generalized) coordinates. 
Note that this is possible at imposed t, but cannot hold 
at fixed 7. Assuming then pair interactions and using 
Eq. one can reformulate Eq. ([55]) within a few lines. 
Without further approximation this yields 



MfI^ 



Of 



x.l 



dxi 



(54) 



where the sum runs now over all interactions / between 
particles i < j. xi and yi refer to components of the 
distance vector of the interaction I and fx.i to the x- 
component of the central force f{ri) = —u'{ri) between 
a particle pair at a distance r; . Note that in Eq. ((54)) the 
stress fluctuations stemming from different interactions 
are decoupled, i.e. MfI^. characterizes the self- or two- 
point contributions of directly interacting beads. Since 
fxi^) = —u'{r)x/r and, hence. 



. dfxjr) 
dx 



{r'^u"{r) — r'^u'(r)) n^n^ 
ru'{r) n^Uy, (55) 



this is identical to the affine shear elasticity ma derived 
at the end of Sec. Ill Bl We have thus confirmed Eq. ([2T|) . 
Interestingly, the argument can be turned around: 

(1) The stress fluctuation contribution term /xp]^ = ma 
is obtained by the simple derivation given in this 
paragraph. 

(2) Using the general Legendre transformation, 
Eq. (f20| . this implies Grr = Ma — Mf for the 
NV7T- and the NP7T-ensemble. 

As already stressed, our thermodynamic derivation of 
Eq. ^ is rather general, not using in particular a well- 
defined displacement field for the particle positions. 
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D. Impulsive corrections for truncated potentials 

Truncation. It is common practice in computational 
condensed matter physics [6, 9] to truncate a pair inter- 
action potential u{r), with r being the distance between 
two particles i and j, at a conveniently chosen cutoff Tc. 
This allows to reduce the number of interactions to be 
computed and energy or force calculations become 0{N)- 
processes. However, the truncation also introduces tech- 
nical difficulties, e.g., instabilities in the numerical so- 
lution of differential equations as investigated especially 
for the MD method [5^ . Without restricting much in 
practice the generality of our results, we assume below 
that 

• the pair potential u(r) is short-ranged, i.e. that it 
decays within a few particle diameters, 

• it scales as u{r) = u(s) with the reduced dimen- 
sionless distance s — r/ai where tr/ characterizes 
the length scale of the interaction I and 

• the same reduced cutoff Sc = r^jai is set for all 
interactions I. 

For instance, for monodisperse particles with constant 
diameter tr, as for the standard Lennard-Jones (LJ) po- 
tential 0, 

ziLj(.)-4e(^-L- J_^, (56) 

the scaling variable becomes s — r/a and the reduced 
cutoff Sc = fc/cr. The effect of introducing Sc is to replace 
u{s) by the truncated potential 

Ut{s)^u{s)H{s,~ s) (57) 

with H{s) being the Heaviside function [4^. Even if 
Eq. (|57p is taken by definition as the new system Hamil- 
tonian, it is well known that impulsive corrections at the 
cutoff have to be taken into account in general for the ex- 
cess pressure Pox and other moments of the first deriva- 
tives of the potential [Q] . This is seen by considering the 
derivative of the truncated potential 

u[{s) = u'{s)H{s^ ~s)~ u{s)5{s - Sc). (58) 

Shifting of truncated potential. These corrections 
with respect to first derivatives can be avoided by con- 
sidering a properly shifted potential [9| 

Us(s) = (u(s) - u(sc)) i/(sc - s) (59) 

since u'^{s) = u' [s)H{sc — s). With this choice no impul- 
sive corrections arise for moments of the instantaneous 
shear stress f and of other components of the excess 
stress tensor Eq. (^5]) . Specifically, if the poten- 

tial is shifted, all impulsive corrections are avoided for 
the shear stress fluctuations /ip, Eq. 



Correction to the Born coefficient. As shown in 
Ref. [3^, the standard shifting of a truncated potential 
is, however, insufficient in general for properties involving 
second (and higher) derivatives of the potential. This is 
particulary the case for the Born coefhcient ^-q , Eq. (|47l) , 
which is required to compute the affine shear elasticity 
fiA, Eq- O- The second derivative of the truncated and 
shifted potential being 

<(s) = u"(s)i/(sc - s) - u'{s)S{s - Sc) (60) 

this implies that the Born coefhcient reads /^b = — 
A/iB with fiB being the bare Born coefficient and A^b 
the impulsive correction which needs to be taken into 
account. The latter correction may be readily computed 
numerically from the configuration ensemble using [3^ 

A^B = lini /ib(s) with 

being a weighted radial pair distribution function which 
is related to the standard radial pair distribution function 
5(r) 

Mixtures and polydisperse systems. Below we shall 
consider model systems for mixtures and polydisperse 
systems where ai may differ for each interaction I. For 
such mixed potentials m(s), Mt(s) and Us{s) and their 
derivatives take in principal an explicit index I, i.e. one 
should write ui{s), Mt,i(s), Ws,i(s) and so on. (This is of- 
ten not stated explicitly to keep a concise notation.) For 
example one might wish to consider the generalization of 
the monodisperse LJ potential, Eq. (I56p . to a mixture or 
polydisperse system with 

ui{s) = Aei (s"^2 - s"^) with s = r/ai (62) 

where ei and ai are fixed for each interaction /. In prac- 
tice, each particle i may be characterized by an energy 
scale Ei and a "diameter" Di. The interaction param- 
eters ei{Ei,Ej) and ai{Di, Dj) arc then given in terms 
of specified functions of these properties The ex- 

tensively studied Kob-Andersen (KA) model for binary 
mixtures of beads of type A and B [ll], is a particu- 
lar case of Eq. with fixed interaction ranges ctaa, 
ctbb and uab and energy parameters eaa, ^bb and cab 
characterizing, respectively, AA-, BB- and AB-contacts. 
The expression Eq. (j6ip remains valid for such explicitly 
Z-dependent potentials. 

Shear modulus and compression modulus. Since 
Grr — ^J■B — Pcx — f-F , Eq. (pT|) implies in turn 

Grr = Gtt — A/iB (63) 

with Gtt being the uncorrected bare shear modulus. As 
we shall see in Sec. IIV A[ the correction A/iB is of impor- 
tance for the precise determination of Gtt close to the 
glass transition where the shear modulus must vanish. 
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An impulsive correction has also to be taken into account 
for the compression modulus K = Kpp computed from 
the normal pressure fluctuations in a constant volume en- 
semble as discussed in Appendix |3 Using t he sy mmetry 
of isotropic systems, Eq. (jA24[) . one verifies |39| 



0.96.Lr 



pLJ N=10 , P=2,T=0 



VP — ^^pp 



2 + d 



(64) 



with Kpp being the uncorrected compression modulus. 



III. ALGORITHMIC AND TECHNICAL ISSUES 

Introduction. The computational data presented in 
this work have been obtained using two extremely well 
studied models of colloidal glass-formers which are de- 
scribed in detail elsewhere [l^ H, E, [13 • We present 
first both models and the molecular dynamics (MD) and 
Monte Carlo (MC) [1, i, ^ methods used to compute 
them in the different ensembles (NV7T, NVrT, NP7T 
and NPtT) compared in this study. We comment then on 
the quench protocol and locate the glass transition tem- 
perature Tg for both models via dilatometry, as shown 
in Fig. [21 This is needed as a prerequisite to our study 
of the elastic behavior in Sec. IIVI As sketched in Fig. [1] 
periodic boundary conditions are applied in all spatial di- 
rections with Lx = Ly = . . . for the linear dimensions of 
the d-dimensional simulation box. Time scales are mea- 
sured in units of the LJ time tlj — (mcr^/e)^/^ IQ] for our 
MD simulations (to being the particle mass, a the refer- 
ence length scale and e the LJ energy scale) and in units 
of Monte Carlo Steps (MCS) for our MC simulations. LJ 
units are used throughout this work (e = cr = to = 1) 
and Boltzmann's constant fee is also set to unity. 

Koh- Andersen model. The Kob-Andersen (KA) 
model [5^ for binary mixtures of LJ beads in d = 3 di- 
mensions has been investigated by means of MD simula- 
tion taking advantage of the public domain LAMMPS 
implementation [SS*]. We use N = n^ + n^ = 6912 beads 
per simulation box and molar fractions Cq = np^/n = 0.8 
and Cf, = n-Q/n — 0.2 for both types of beads A and B. 
Following Ref. ^] we set caa — l-Ou, ctbb ~ 0.88ct and 
Cab = 0.8(7 for the interaction range and eAA = l.Oe, 
Ebb = 0.5e and eAB = l-5e for the LJ energy scales. 
Only data for the usual (reduced) cutoff Sc = r^joi = 2.5 
are presented. For the NP7T-ensemble we use the 
Nose-Hoover barostat prov ided by the LAMMPS code 
("fix npt command") [33| which is used to impose a 
constant pressure P = 1 for all temperatures. Starting 
with an equilibrated configuration at Tg = 0.7 well 
above the glass transition temperature Tg — 0.41 (as 
confirmed in the inset of Fig. [2|), the system is slowly 
quenched in the NP7T-ensemble. The imposed mean 
temperature T{t) varies linearly as T{t) = Ti — Rt with 
R — 5 X 10~^ being the constant quench rate. After a 
tempering time icq — we compute over a sampling 



" R=10 

° R=10"' 

O R=10"° 

^ R=10 "^ 

O R=10"° 

X R=10"'forNPiT 

^ final best value 




time t — 5 ■ 10 various properties for the isobaric 
system at fixed temperature such as the moduli Grr 



FIG. 2: Density p{T) and determination of glass transition 
temperature Tg. Inset: Plotting log(l/p) vs. T for the KA 
model as suggested in Ref. Q reveals a linear low-T (solid 
line) and a linear high-T regime which confirms Tg « 0.41 [s^] . 
Main panel: Number density p for pLJ systems at constant 
normal pressure P = 2 as a function of temperature T for 
different quench rates R in units of MCS. The crosses refer to a 
quench in the NPrT-ensemble where the box shape is allowed 
to fluctuate at r = 0. All other data refer to the NP7T- 
ensemble with 7 = 0. The solid line and the dashed line are 
linear fits to the best density estimates (stars) obtained for, 
respectively, low and high temperatures. By matching both 
lines one determines Tg « 0.26. 



represented by the filled spheres in Fig. [TOl Fixing only 
then the volume (NV7T-ensemble) the system is again 
tempered (toq = 10^) and various properties as the 
ones recorded in Table |T] are computed using a sampling 
time t = 5 • 10^. Note that in the NV7T-ensemble 
the equations of motion are integrated using a velocity 
Verlet algorithm [9J with a time step St = 0.01 and 
systems are kept at the imposed temperature T using a 
Langevin thermostat ;9] with friction coefficient 7 = 0.5. 
By redoing for several temperatures around and below 
Tg the tempering and the production runs we have 
checked that the presented results remain unchanged 
(within numerical accuracy) and that ageing effects 
are absolutely negligible. Averages are taken over four 
independent configurations. As one may expect from 
the discussion in Sec. Ill C[ very similar results have 
been found for both ensembles for simple averages, for 
instance for the affine shear elasticity /ia, and for the 
shear modulus Grr, Eq- ©■ Unfortunately, we are 
yet unable to present data for the KA model data at 
imposed shear stress r (NVrT-, NPrT-ensembles). 

Polydisperse LJ beads. A systematic comparison of 
constant-7 and constant-T ensembles has been per- 
formed, however, by MC simulation of a specific case 
of the generalized LJ potential, Eq. ([5^ . where all inter- 
action energies are identical, e; = e, and the interaction 
range is set by the Lorentz rule cr/ — (Di + Dj)/2 [sij 
with Di and Dj being the diameters of the interacting 
particles. Only the strictly two-dimensional (2D) case 
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TABLE I: Some properties for the KA model at normal pres- 
sure P = 1 and shear stress r = 0: the mean density p, the 
mean energy per volume ep, the impulsive correction term 
Ap,B (Sec. Ill pp . the (corrected) Born-Lame coefficient pB, 
the afhne shear elasticity pA ~ Pb — Pcx, the shear modulus 
G-TT obtained in the NV7T-ensemble, the hypervirial tjb and 
the affine dilatation elasticity rjA = 2Pid + t;b + Pcx discussed 
in Appendix [X] and the compression modulus Kpp obtained 
for the NV7T-ensemble. The affine elasticities pA and rjA are 
the natural scales for, respectively, the shear modulus G and 
the compression modulus K. 

{d — 2) is considered. Following Ref. the bead diam- 
eters of this polydisperse LJ (pLJ) model are uniformly 
distributed between O.Scr and 1.2a. For the examples re- 
ported in Sec.|TV]we have used N — 10000 beads per box. 
We only present data for a reduced cutoff Sc = 2.0so with 
So — 2^/^ being the the minimum of the polydisperse LJ 
potential, Eq. (|62p . Various properties obtained for the 
pL J model are summarized in Table HIl 

Local MC moves for the pLJ model. For all ensem- 
bles considered here local MC moves are used (albeit not 
exclusively) where a particle is chosen randomly and a 
displacement Sr, uniformly distributed over a disk of ra- 
dius (Sr,„ax, from the current position r of the particle 
is attempted. The corresponding energy change SE is 
calculated and the move is accepted using the standard 
Metropolis criterion psj . The maxium displacement dis- 
tance <5rmax is chosen such that the acceptance rate A 
remains reasonable, i.e. 0.1 < A < 0.5 As may be 
seen from the inset in Fig. [3] where the acceptance rate 
A{T) is shown as a function of temperature T, we find 
that temperatures below T = 0.1 are best sampled us- 
ing (5r„iax ~ 0.01, whereas Jrmax « 0.1 is a good choice 
for the interesting temperature regime between T = 0.2 
and T = 0.3 around the glass transition temperature 
Tg = 0.26. (See the main panel of Fig. [5] for the de- 
termination of Tg.) Various data sets are presented in 
Sec. llVl for a fixed value (5rmax = 0.1 which allows a rea- 
sonable comparison of the sampling time dependence for 
temperatures around the glass transition. This value is 
not necessarily the best choice for tempering the system 
and for computing static equilibrium properties at the 
given temperature, especially below T « 0.2. 
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TABLE 11: Some properties of the pLJ model at P — 2 and 
T = 0: the mean density p, the mean energy per volume ep, 
the impulsive correction term ApB, the affine shear elasticity 
PA, the shear moduli G^t and G^-y obtained in the NPrT- 
ensemble and Gtt obtained in the NV7T-ensemble, the affine 
dilatation elasticity t/a, the compression modulus K^p ob- 
tained in the NPrT-ensemble and the Rowlinson formula Kpp 
for the NV7T-ensemble. All data have been obtained for a 
sampling time t = 10^ MCS. Please note that the values given 
correspond to the best values obtained with the most appro- 
priate technique and may differ thus slightly from the data 
represented in the figures for a specific protocol. 



Collective plane-wave MC moves. Local MC jumps in 
the NV7T-ensemble become obviously inefficient for re- 
laxing large scale structural properties f25l| for the lowest 
temperatures computed for the pLJ model (T < 0.1) as 
may be seen, e.g., from the finite, albeit very small shear 
stresses r which happen to occur naturally in some of the 
quenched configurations (especially if the quench rate R 
is too large) and can o nly with difficulty be relaxed us- 
ing local jumps alone [53|. We have thus crosschecked 
and improved the values obtained using only local MC 
moves by adding global MC move attempts for all parti- 
cles with longitudinal and (more importantly) transverse 
plane waves commensurate with the simulation box. The 
random amplitudes and phases of the waves are assumed 
to be uniformly distributed. Note that the maximum 
amplitude aniax(9) of each wave is chosen inversely pro- 
portional to the length g = |g| of the wavevector in agree- 
ment with continuum theory. This allows to keep the 
same acceptance rate for each global mode if g ^ 1. As 
one expects, the plane- wave displacement attempts be- 
comes inappropriate for large q where continuum elastic- 
ity breaks down [l3, [13, 13 Si,nd, hence, the acceptance 
rate A{q) gets too large to be efficient [sJl. Note that 
the maximum amplitude ainax('?) has to be chosen much 
smaller for the longitudinal waves as for the transverse 
ones, i.e. since the systems are essentially incompress- 
ible, the transverse modes naturally are more effective for 
quenching and equilibrating the systems. These global 
plane-wave moves have been added for the tempering of 
all presented configurations below T = 0.3 with tcq = 10^ 
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FIG. 3: Energy per particle e (main panel) and acceptance 
rate A (inset) obtained for pLJ systems quenched at constant 
pressure P = 2 at a given quench rate R as indicated. The 
crosses refer to a quench at imposed t = 0. The solid line and 
the dashed line are linear fits to the energy for, respectively, 
low and high temperatures confirming the value Tg 0.26 
for the pLJ model. Only local MC moves with one fixed 
(non-adaptive) maximum particle displacement 5r-niax for all 
temperatures T have been used for the MC data shown in 
the main panel. The acceptance rate decreases thus with 
decreasing T as shown in the inset for (5rmax ~ 0.01 (upper 
data) and Sr^iax = 0.1. 



and for production runs over sampling times t = 10 for 
static properties for T < 0.2. Details concerning the col- 
lective plane-wave MC moves must be given elsewhere. 

MC sampling of different ensembles. In this paper 
we shall be concerned more with the consequences of 
afBne displacements associated with collective MC moves 
changing the overall box volume V and shear strain 7. 
Using again a Metropolis criterion for accepting a sug- 
gested change of V or 7, these pure dilatational and pure 
shear strain fluctuations are used to impose a mean nor- 
mal pressure P — 2 and a mean shear stress t — 0. For a 
shear strain altering MC move one first chooses randomly 
a strain fluctuation ^7 with \S^\ < 57max- Assuming an 
imposed shear stress r the suggested ^7 is accepted if 

^ < exp {-PSE + PSjt) (65) 

with ^ being a uniformly distributed random variable 
with < ^ < 1. The energy change 5E associated with 
the afiine strain may be calculated using Eq. (I55t . Since 
T = is supposed in this study, the last term in the ex- 
ponential drops out. (As noted in Sec. IIIBl the volume 
remains unchanged in a pure shear strain and, hence, no 
translational entropy contribution appears.) The max- 
imum shear strain step Sj^i^^ is fixed such that at the 
given temperature the acceptance rate remains reason- 
able dl]. As described in more detail in Ref. Q, one 
similarly chooses a relative volume change 5V/V with V 
being the current volume [55| . The volume altering move 



is accepted at an imposed pressure P if 

C < exp i-l36E + N log(l + SV/V) - 135V P) (66) 

where 5E may be computed using Eq. (|A16p . The log- 
arithmic contribution corresponds to the change of the 
translational entropy. Using Eq. ((65)) after each MCS 
for local MC moves (and global plane-wave MC moves if 
these are added) realizes the NVrT-ensemble shown in 
Fig. IHc), using Eq. the NP7T-ensemble shown in 
Fig. [ijd) and using in turn both strain fluctuations the 
NPrT-ensemble in Fig. (He). We remind that NV7T- and 
NP7T-ensembles are used to determine Grr and NVrT- 
and NPrT-ensembles to obtain G^r and G^^. As shown 
in Fig. [31 we quench the configurations starting from 
Ti = 0.5 in the NP7T- and the NPrT-ensemble using 
a constant quench rate R. Note that the smallest quench 
rate R = 10~* for one configuration in the NP7T-shown 
did require alone a run over six months using one core 
of an Intel Xcon E5410 processor. Interestingly, similar 
results are obtained with the NPrT-ensemble only using 
a rate R = 10~^. (Due to the additional computation 
for the shear strain fluctuations this is, however, only 
a factor flve faster.) The configurations created by the 
NP7T-quenches are used (after tempering) for the sam- 
pling in the NV7T- and NP7T-ensembles, the configura- 
tions obtained by the NPrT-quenches for the sampling 
in the NVrT- and NPrT-ensembles. Only two indepen- 
dent configurations have been sampled following the full 
protocol for each temperature and each of the four en- 
sembles. Instead of increasing further the number of con- 
figurations (which will be done in the future) we have 
focused in the present preliminary study on long temper- 
ing times and production runs over t — 10^ MCS which 
have been redone several times for a few temperatures 
(summing up to total runs of « 5 • 10^ MCS for T = 0.2 
and T = 0.25) to check for equilibration problems and 
to verify that ageing effects can be ignored. We firmly 
believe that the breaking of the time translational invari- 
ance is not a relevant issue for the strong sampling time 
dependence reported below for the shear moduli. 

Glass transition temperature. The inset of Fig. [2] 
presents for the KA model how the glass transition tem- 
perature Tg may be obtained using standard dilatometry 
where the rescaled density fog(l/p) is traced as a func- 
tion of temperature T as in Ref. [7| for a similar polymer 
glass problem. This reveals a linear low-T (solid line) and 
a linear high-T regime. The intercept of both lines con- 
firms the well-known value Tg « 0.41 from the literature 
[5^ . The main panel presents the unsealed number den- 
sity p for pLJ systems at constant normal pressure P = 2 
as a function of temperature T for different quench rates 
R in MCS. Only local MC moves with a fixed maximum 
particle displacement Sr-a^ax = 0.01 for all temperatures 
T are used for the examples given. The crosses refer to 
a quench in the NPrT-ensemble where the box shape is 
allowed to fluctuate at r = 0. All other data refer to 
the NP7T-ensemble with 7 = 0. The solid line and the 
dashed line are linear fits to the best density estimates 
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FIG. 4: Rescaled shear modulus G{t)/fiA vs. sampling time 
t in the high-T liquid limit of the pLJ model comparing var- 
ious ensembles and observables. The filled squares show the 
bare shear modulus Gtt obtained using Eq. ((3]) for the NV7T- 
ensemble without the impulsive correction A/ib ~ 0.2 (hori- 
zontal line). The corrected modulus Gtt ~ Gtt — A/iB for 
the NV7T-ensemble (open squares) and the NP7T-ensemble 
(open diamonds) vanishes as it should. Note that Gtt de- 
creases slightly faster than G-yi- ~ G-y^. As shown by the 
dashed line, all operational definitions of G decay inversely 
with t. The dash-dotted line indicates the expected power- 
law exponent —1/2 for the correlation coefficient c^t (stars). 



indicated by the stars obtained for, respectively, low and 
high temperatures. By matching both lines one deter- 
mines Tg 0.26. A similar value is given if log(l/p) is 
plotted as a function of T. The main panel of Fig. [3] 
shows the interaction energy per particle e of the pLJ 
model as a function of temperature for the same quench 
protocols as in Fig. [H The solid line and the dashed 
line are linear fits to the energy for, respectively, low and 
high temperatures. This confirms Tg 0.26 for the pLJ 
model 15611. 



IV. COMPUTATIONAL RESULTS 

A. High-temperature liquid limit 

We begin our discussion by focusing on the high tem- 
perature liquid limit where all reasonable operational 
definitions of the shear modulus G must vanish. The 
results presented in Fig. 0] have been obtained for the 
pLJ model at mean pressure P — 2, mean shear stress 
T = and mean temperature T — 1, i.e. far above the 
glass transition temperature Tg w 0.26. Only local MC 
monomer displacements with a maximum jump displace- 
ment (5riiiax = 0.1 are reported here. Instantaneous prop- 
erties relevant for the moments are written down every 
10 MCS and averaged using standard gliding averages Q, 
i.e. we compute mean values and fluctuations for a given 



time interval [ioi = + ^] and average over all possible 
intervals of length t. The horizontal axis indicates the 
latter interval length t in units of MCS. The vertical axis 
is made dimensionless by rescaling the moduli with the 
(corrected) affine shear elasticity fj,A = ^J■B — fox ~ 17 
from Table ini (Being a simple average, fiA is found to be 
identical for all computed ensembles.) 

In agreement with Eq. ([^5]) . the ratio G{t)/iiA for all 
operational definitions drops rapidly below unity. The 
dashed line indicates the asymptotic power-law slope — 1 
[2l] . [soj . This can be understood by noting that at im- 
posed T the shear strain freely diffuses in the liquid limit, 
i.e. (^7^) ~ t, which implies that G « G^yr ~ 1/t. For 
both NVtT- and NPrT-ensembles we compare our key 
definition G-yr, Eq. to the observable Eq. In 
both cases wc confirm G-yr (t) ~ G-y-y (t) for larger times as 
suggested by the thermod yna mic relation Eq. ([5]) which 
has been directly checked |38| . 

We have also computed the dimensionless correlation 
coefficient Cyrit), Eq. ([7]), which characterize the corre- 
lation of the measured instantaneous shear strains 7 and 
shear stresses f. As indicated for the NPrT-ensemble 
(stars), Cyr vanishes rapidly with t, i.e. 7 and f are decor- 
related. The dash-dotted line indicates the power-law ex- 
ponent — 1/2 expected from Eq. and Gyr{t) ~ 1/t. 
Note that by plotting c-yrit) vs. G-yr{t)/nA, Eq. ([22]) 
has been directly verified for both NVrT- and NPrT- 
ensembles (not shown). 

The stress fluctuation formula, Eq. ([3]), is represented 
in Fig. m by squares and diamonds for, respectively, 
NV7T- and NP7T-ensembles. The bare shear modulus 
Gtt without the impulsive correction A/xb for the Born 
coefficient /ie discussed in Sec. Ill Dl is given by small filled 
symbols. The impulsive correction A/ib ~ 0.2 (Table HI]) 
computed independently from the weighted radial pair 
correlation function, Eq. (j61|) . is shown by a horizontal 
line. The corrected moduli Gtt — Gtt — A/ib vanish in- 
deed as expected. See Ref. [39| for a systematic numerical 
characterization of A/iB as a function of the reduced cut- 
off Sc. As may also be seen there, a truncation correction 
is also necessary for the KA model. The A/ib for dif- 
ferent temperatures may be found in Table [T] for the KA 
model and in Table |TT] for the pLJ model. We assume 
from now on that the impulsive corrections are properly 
taken into account without stating this technical point 
explicitly. 



B. Dynamical matrix harmonic network 

Introduction. We turn now to the opposite low-T 
limit focusing specifically on the pLJ model at T = 0.001. 
Obviously, such deeply quenched colloidal glasses must 
behave as amorphous solids with a constant shear modu- 
lus for large times (albeit not infinite times) as sketched 
in Fig. [TJb) by the top curve. Before presenting in 
Sec. IIVCI the elastic properties of these glasses, we dis- 
cuss first conceptually simpler substitute systems formed 
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by permanent spi'mg networks. We do this for illustration 
purposes since questions related to ergodicity and ageing 
are by construction irrelevant and the thermodynamical 
relations reminded in Sec. HIl should hold rigorously. 

Network Hamiltonian. Assuming ideal harmonic 
springs the interaction energy reads 



(67) 



with Ki being the spring constant and Ri the reference 
length of spring I. The sum runs over all springs / be- 
tween topologically connected vertices i and j of the net- 
work. We assume here a strongly and homogeneously 
connected network where every vertex i is in contact with 
many neighbors j. Such a network may be constructed 
from a pLJ configuration at T = 0.001 keeping the parti- 
cle positions for the positions of the vertices and replacing 
each LJ interaction I by a spring of spring constant Ki 
and reference length Ri. The simplest choice we have 
investigated is to set Ki = 1 and Ri = ri which corre- 
sponds to a well-defined ground state of energy E = Q 
and pressure P = at T = (not shown). We discuss 
here instead a slightly more realistic choice of Ki and 
Ri which corresponds to the same "dynamical matrix" 
as the pLJ configuration at T and P — 2, i.e. the 
same second derivative of the total interaction potential 
with respect to the particle positions [l3,[i3]- This choice 
imposes the setting 



Ri 



u'{r) 



n 



u"{r) 



(68) 



with u{r) = u{s) being the interaction potential for re- 
duced distances s/ = ri/cri < Sc- (Impulsive corrections 
at si = Sc are neglected here.) 

Sampling. Taking advantage of the fact that all in- 
teracting vertices are known, these networks can be com- 
puted using global MC moves with longitudinal and 
transverse planar waves as described in Sec. lIIIl No local 
MC jumps for individual vertices have been added. Note 
that although some Ki and Ri may even be negative, this 
does not affect the global or local stability of the network. 
Presumably due to the fact that our spring constants are 
rather large, we have not observed either any buckling in- 
stability if volume fluctuations are allowed. For networks 
at the same pressure P, shear stress r and temperature T 
as the original pLJ configuration, we obtain, not surpris- 
ingly, the same afhne shear elasticity /xa = 33.9 as for the 
reference. The same applies for other simple averages. 

Shear modulus. Figure [5] presents various measure- 
ments of the shear modulus G for different ensembles as 
a function of sampling time t. As in Sec. IIV Al we use 
gliding averages and the vertical axis has been rescaled 
by the afhne shear elasticity ^a- As can be seen, we 
obtain in the long-time limit G//iA ~ 0.67 (bold line) 
irrespective of whether the shear modulus is determined 
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FIG. 5: Rescaled shear modulus G{t)/^,A vs. sampling time 
t for systems of permanent linear springs corresponding to 
the dynamical matrix of pLJ systems at P = 2, r = and 
T — 0.001. As indicated by the bold line, we obtain G/ha ~ 
0.67 in the long-time limit. In agreement with Eq. (|21|) Gtt 
vanishes with time if it is computed in an ensemble where 7 
fluctuates as shown for the NVrT-ensemble (filled squares). 



from G-yr or G-y^ in the NVtT- or NPrT-ensembles or 
using Grr in the NV7T- or NP7T-ensembles at fixed 
shear strain 7. Note that even in these cases the shear 
modulus decreases first with time albeit the network is 
perfectly at thermodynamic equilibrium and no ageing 
occurs by construction. As shown by the dash-dotted 
line, G(<)//iA ^ 1 for short times, i.e. in this limit the 
response is affine. As indicated by the stars, we have also 
computed the correlation coefficient Cjr(t) as a function 



of time. From c-, 



1 for small t, it decreases somewhat 
0.82 for large times in agree- 



becoming Cjr w VO.67 
ment with Eq. (|22|) . If Gtt is computed in the "wrong" 
NVtT- or NPrT-ensembles, it is seen to vanish rapidly 
with time as shown by the filled squares. This decay is 
of course expected from Eq. ((?T|l as discussed in Sec. Ill Bl 
and Sec. IXTOl 

Shear-stress fluctuations. To emphasize the latter 
point we show in Fig. [5] the stress fluctuations Hvit) for 
different ensembles. This is also done to have a closer 
look at the time dependence of various contributions to 
G'T-T-(t). As before, the vertical axis has been rescaled 
with /iA = 33.9. We also indicate the rescaled values of 
/iB(0 (thin top line) and HAit) = ^Bit) — Pex(i) (dash- 
dotted line) computed by a gliding average over time 
intervals [^0,^1 = to + t]. Being simple averages both 
properties do not depend on the ensemble probed (not 
shown) and become virtually immediately i-independent 
as can be seen from the indicated horizontal lines. This 
is quite different for the different shear stress fluctuations 
/ip presented which increase monotonously from /xp ~ 
for short times to their plateau value for long times. In 
agreement with Eq. (pTt we find that fJ.F{t)\^ — 7> ^a 
for NVtT- and NPrT-ensembles and as predicted by 
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FIG. 6: Shear stress fluctuations /iF(f)//^A with /ia = 33.9 
for the same systems as in Fig. [5] We observe /ip(t)/^A — > 
1 (dash-dotted hne) for imposed r, while ^F{t)/jj,A 1 — 
G//XA ~ 0.32 (dashed hne) for imposed 7. Note the broad 
time-dependence of fJ,F{t) in the latter cases. 
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FIG. 7: G^T{t)/iJ,A, Gj^{t)/jj,A and GTTit)/fJ,A for pLJ beads 
in the low-temperature limit for P — 2, t = and T = 
0.001 where /ia ~ 33.9. Qualitatively, the data compares well 
with the various shear moduli presented in Fig. [S] The filled 
triangles show Gtt for the NPrT-ensemble which is seen to 
vanish rapidly with sampling time t. The dashed line indicates 
a 1/t-decay. 



Eq. dini) we confirm /iF(i)|^ ha - G for NV7T- and 
NP7T-ensembles. The asymptotic limit is reached after 
about lO'' MCS in the first case for imposed t, but only 
after 10^ for imposed 7. This is due to the fact that ^fI,- 
corresponds to the fluctuations of the self-contribution 
of individual springs as discussed in Sec. Ill Cl while /ip]^ 
also contains contributions from stress correlations of dif- 
ferent springs. It should be rewarding to analyze in the 
future finite system-size effects and the time dependence 
of three-dimensional (3D) systems using the dynamical 
matrix associated to the KA model in the T — > limit. 



C. Low-temperature glass limit 

We return now to the colloidal glasses where the topol- 
ogy of the interactions is, of course, not permanently 
fixed (at least not at finite temperatures) and the shear 
stresses may thus fiuctuate more strongly. As an exam- 
ple, we present in Fig.[7]the shear modulus G{t) as a func- 
tion of sampling time t for pLJ beads at mean pressure 
P = 2, mean shear stress t = and mean temperature 
T — 0.001, i.e. the same conditions as in the previous 
subsection. All simple averages, such as the affine shear 
elasticity /ia — 33.9, are the same for all ensembles albeit 
different quench protocols have been used. For all exam- 
ples given we use local MC jumps with frmax = 0.01 
which allows to keep a small, but still reasonable ac- 
ceptance rate A « 0.161. We have also performed runs 
with additional global MC moves using longitudinal and 
transverse plane waves. This yields similar data which is 
shifted horizontally to the left (not shown). 

As one would expect, the shear moduli presented in 
Fig. [7] are similar to the ones shown in Fig. [5] for the per- 
manent spring network. The filled triangles show Grr 
computed in the "wrong" NPrT-ensemble. As expected 
from Eq. ([^0)1 these values vanish rapidly with time. If on 
the other hand G-yr, Gjj and Grr are computed in their 
natural ensembles, all these measures of G are similar and 
approach, as one expects for such a low temperature, the 
same finite asymptotic value G/ ~ 0.46 indicated by 
the bold horizontal line. The stars indicate the squared 
correlation coefficient c^T-(i) which is seen to collapse per- 
fectly on the rescaled modulus G^r{t) I y^K which shows 
that Eq. even holds for small times before the ther- 
modynamic plateau has been reached. Similar reduced 
shear moduli Grr/ Ma have been also obtained for the KA 
model in the low-T limit as may be seen from Table [D 
These results confirm that about half of the affine shear 
elasticity is released by the shear stress fluctuations in 
agreement with Refs. [13, [3 ■ Interestingly, this value is 
slightly smaller as the corresponding value for the per- 
manent dynamical matrix model presented above. Since 
some, albeit very small, rearrangement of the interactions 
must occur for the pLJ particles at finite temperature, 
this is a reasonable finding. 

The rescaled shear stress fluctuations iJL-p{t) / l^-A for 
the same ensembles are presented in Fig. [S] While the 
two-point correlations fiB{t)/lJ-A (thin horizontal line) 
and /iA(i)//^A (dash-dotted horizontal line) do not de- 
pendent on time, one observes again that Hp(t)/fj,A in- 
creases monotonously from zero to its long-time asymp- 
totic plateau value. In agreement with Eq. (j2ip we find 
fJ-F{t)\^ ^J-A for the NVrT and the NPrT ensembles, 
while iji.F{t)\^ HA — G (dashed horizontal line) for the 
NV7T and the NP7T ensembles. (Similar behavior has 
been obtained for the KA model in the NV7T-ensemble 
for temperatures T < 0.1.) As may be better seen for 
the data set obtained for T = 0.010 using the NV7T- 
ensemble (small filled squares) where fiA ~ 33.7, the 
stress fiuctuations /iF(t)L do not rigorously become con- 
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FIG. 8: Reduced shear stress fluctuation /ip(f)//iA for pLJ 
beads using the same representation as in Fig. [S] The small 
filled squares refer to data obtained for T — 0.010 using the 
NV7T-ensemble, all other data to T = 0.001. The reduced 
Born coefficient fiB{t)/l^A (thin line) and the affine shear elas- 
ticity fiA{t)/^J.A (dash-dotted line) are essentially time inde- 
pendent while /ip(t)/^A approaches its large-t limit (dashed 
line) from below. 



slant, but increase extremely slowly on the logarithmic 
time scales presented. There is, hence, even at low tem- 
peratures always some sampling time dependence if the 
interactions are not permanently fixed as for the topologi- 
cally fixed networks discussed above. However, this effect 
becomes only sizeable above T = 0.2 for the KA model 
and above T = 0.1 for the pLJ model. We shall now 
attempt to characterize this temperature dependence. 



D. Scaling with temperature 

Stress fluctuations. We turn now to the characteriza- 
tion of the temperature dependence of the shear modu- 
lus G and various related properties. We focus first on 
the best values obtained for the longest simulation runs 
available. Figure IHl presents the stress fluctuations ob- 
tained for the KA model using the NV7T-ensemble (open 
spheres) and for the pLJ model using NV7T- and the 
NPrT-ensembles. The affine shear elasticity ha obtained 
for each system is represented by small filled symbols. As 
shown for both models by a bold solid line for the low-T 
and by a dashed line for the high-T regime, IjLa{T) de- 
cays roughly linearly with T . Interestingly, the slopes 
match precisely at T = Tg for both models. As may be 
seen from the large triangles for the NPrT-ensemble, we 
confirm for all temperatures that /ip]^ ~ /iA as expected 
from Eq. pT|) . For systems without box shape fluctua- 



tions, it is seen that /ipl-^ becomes non-monotonous with 
a clear maximum at T w Tg. In the liquid regime for 
T > Tg where the boundary conditions do not matter, 



we obtain ji-p] 



Mf| 



fiA, i-e. all shear stress fluctu- 



ations decay with temperature. Below Tg the boundary 
constraint (7 — 0) becomes more and more relevant with 
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FIG. 9: Stress fluctuation hf (open symbols) and affine shear 
elasticity (small flUed symbols) vs. temperature T for 
both models. The indicated values have been obtained from 
the longest simulation runs available for a given temperature. 
fiA decays roughly linearly with T as shown by bold solid 
lines for the low-T and by dashed line for the high-T regime. 
For constant-r systems and for all systems with T > Tg we 
confirm that /ip ~ /ia, while fipl^ (open spheres and squares) 
is seen to be non-monotonous with a maximum at T^. 
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FIG. 10: Unsealed shear modulus G as a function of tem- 
perature T for both models. For the pLJ model the values 
Gtt obtained in the NV7T ensemble (squares) are essentially 
identical to the moduli G^r (crosses) and G^-, (triangles) ob- 
tained for constant t. As shown by the small filled triangles, 
Gtt vanishes in the NPrT-ensemble for all T in agreement 
with Eq. (|21[) . The solid and the dashed line indicate the 
cusp-singularity, Eq. (1691) . for respectively the KA model and 
the pLJ beads 



decreasing T which reduces the shear stress fluctuations. 
According to Eq. (^0]) . /ipj^ must thus decrease with in- 
creasing G as the system is further cooled down. That 
fj-pl^ must necessarily be non-monotonous with a maxi- 
mum at Tg is one of the central results of the presented 
work. 
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FIG. 11: Rescaled shear modulus y — G/ha vs. reduced 
temperature x = T/Tg for the KA model at P = 1 (NV7T- 
ensemble) and the pLJ model at P = 2 (NVrT-, NPtT-, 
NV7T- and NP7T-ensembles). The bold line indicates again 
a continuous cusp-singularity. Inset: Squared correlation co- 
efficient for NVrT- and NPrT-ensembles. 



Shear modulus as a function of temperature. The 
temperature dependence of the (unsealed) shear mod- 
uli G{T) is presented in Fig. [lOl Data obtained by MD 
simulation of the KA model {P ^ 1, Tg 0.41) are in- 
dicated for the NV7T-ensemble (open spheres) and the 
NP7T-ensemble (filled spheres) . All other results refer to 
the best values obtained by MC simulations of the pLJ 
model (P = 2, Tg « 0.26) using both local and global 
moves and different ensembles as indicated. For the pLJ 
model the values Gtt obtained in the NV7T ensemble 
(squares) are seen to be essentially identical for all T to 
the moduli Gjr (crosses) and Gj-y (triangles) obtained 
in ensembles with strain fluctuations. As shown by the 
small filled triangles, Gtt vanishes in the NPtT ensemble 
for all T in agreement with Eq. (I?!]) . (We remind that the 
impulsive correction A/iB has to be taken into account if 
Gtt is used.) Decreasing the temperature further below 
Tg the shear moduli are seen to increase rapidly for both 
models. As shown by the solid and the dashed lines both 
models are well fitted by a cusp-singularity 



G(T) « .91 (1 - r/Tg)i/2 ^ f^^. y ^ J. 



(69) 



with fit constants gi w 23 for the KA model and gi ~ 15 
for the pLJ model where we have set in both cases qn = 0. 
Confirming the MD simulations by Barrat et al. [11| , this 
suggests that the transition is very sharp, albeit contin- 
uous in qua litative agreement with the predictions from 
Ref. j28l. l29l| . If we fit the data close to Tg with an ad- 
ditional off-set go J a slightly negative value is obtained 
which is not compatible with MCT [3l|, [s^]. However, 
admittedly both the number of data points close to Tg 
and their precision (due to the small number of inde- 
pendent configurations) are yet not sufhcient to rule out 
completely a positive, albeit (presumably) small off-set. 
A slightly different representation of the data is given 
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FIG. 12: Unsealed shear modulus G^t as a function of sam- 
pling time t in MCS obtained for the pLJ model in the NPrT- 
ensemble for different temperatures as indicated. The thin 
line indicates the power-law slope —1 for large times in the liq- 
uid regime. A clear shoulder (on the logarithmic scales used) 
can only be seen below T « 0.2 and only below T — 0.01 
one observes a t-independent plateau over two orders of mag- 
nitude. We emphasize that we present here a sampling time 
effect and not a time correlation function. The same sampling 
time effect is observed after additional tempering. 



in Fig. [TT] where the reduced shear modulus y = 
G{T) / iia{T) is plotted as a function of the reduced tem- 
perature X = T/Tg for both models. Using the indepen- 
dently determined glass transition temperature Tg and 
affine shear elasticity /xa(T) as scales to make the axes 
dimensionless, the rescaled data are found to collapse ! 
This is of course a remarkable and rather unexpected re- 
sult considering that two different models in two different 
dimensions have been compared. The bold line indicates 
the cusp-singularity y = 0.45(1 — x)-^^^ with a prefactor 
which is compatible to the ones used in Eq. (jBHl) consider- 
ing the typical values of /xa(T) in both models. Whether 
this striking collapse is just due to some lucky coinci- 
dence or makes manifest a more general universal scaling 
is impossible to decide at present. (Please note that the 
corresponding data for the compression modulus shown 
in Fig. [15] does not scale.) No rescaling of the vertical 
axis is necessary for the squared correlation coefficient 
c^^ shown in the inset for two different ensembles of the 
pLJ model. The same continuous cusp-like decay (bold 
line) is found as for the modulus G. This is again consis- 
tent with the thermodynamic reasoning put forward at 
the end of Sec. MM Eq. (EH). 

Sampling time dependence. A word of caution must 
be placed here. The results presented in Figs. |9][TT] corre- 
spond to the longest simulation runs we have at present 
been able to perform. Obviously, this does not mean that 
they correspond to the limit for asymptotically long sam- 
pling times, if the latter limit exists, or at least an inter- 
mediate plateau of respectable period of time. Focusing 
on the pLJ model in the NPrT-ensemble we attempt in 
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Fig. [12] to give a tentative characterization of the sam- 
phng time dependence. We plot Gjr{t) as a function of 
t for different temperatures as indicated. To make the 
data comparable all simulations for T > 0.1 have been 
performed with local MC jumps of same maximum dis- 
tance Sr jnsix = 0.1. As we have already seen in Fig. |4l 
the shear modulus decays inversely with time in the liq- 
uid regime at high temperatures. This limit is indicated 
by the thin line. With decreasing T the (unsealed) shear 
modulus increases and a shoulder develops. Note that 
even for T = 0.25, i.e. slightly below the glass tran- 
sition temperature Tg = 0.26 estimated via dilatometry 
(Fig. [2]), Gyr{t) decays strongly with time. A more or less 
t-independent shoulder (on the logarithmic scales used) 
can only be seen below T 0.2. We emphasize that we 
present here a sampling time effect expressing the fact 
that the configuration space is more and more explored 
with increasing t and not a time correlation function Q . 
Please note also that this time dependence has nothing 
to do with equilibration problems or ageing effects. Time 
translational invariance is perfectly obeyed in our simula- 
tions as we have checked by rerunning the sampling after 
additional tempering over at least 10^ MCS for aU tem- 
peratures. Similar t-dependencies for different tempera- 
tures T have also been observed for G^r{t) and Grrit) 
for the pL J model and for Grr (t) for the KA model (not 
shown) . 

A different representation of the data for both mod- 
els is shown in Fig. [T3] where the reduced shear modulus 
y = G{T)/ fj,A{T) is plotted against the reduced temper- 
ature X = T/Tg for different sampling times t as indi- 
cated. Similar behavior is found for both models. The 
shear modulus decreases with t and this the more the 
larger the temperature T. For smaller temperatures the 
different data sets appear to approach more readily as 
one expects from the i-independent shoulder for the pLJ 
model in Fig. [TH As a consequence the glass transi- 
tion is seen to become sharper with increasing t. The 
data approach the continuous cusp-singularity indicated 
by the bold line. However, it is not clear from the cur- 
rent data if the data converge to an at least intermedi- 
ately stable behavior. Longer simulation runs are clearly 
necessary for both models to clarify this issue. Similar 
data have also been obtained for the pLJ model from 
the independent determination of G~fr{t) and G^~f{t) in 
the NVrT-ensemble and for Grrit) in the NV7T- and 
NP7T-ensemble . 



V. CONCLUSION 

A. Summary 

The stress fluctuation formalism is a powerful method 
for computing elastic moduli for canonical ensembles 
using simulation boxes of constant volume and shape 
Focusing on the shear modulus G of two well- 
known glass-forming colloidal model systems in d = 3 
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FIG. 13: Reduced shear modulus G(T)/iia{T) as a func- 
tion of reduced temperature x = T/Tg for different sampling 
times t as indicated: (a) GrTiT) for the KA model (NV7T- 
ensemble), (b) G~fT{T) for the pLJ model (NPrT-ensemble). 
For both models the glass transition is seen to become sharper 
with increasing t. The continuous cusp-singularity is indi- 
cated by the bold line. More points in the vicinity of a; ~ 1 
and large sampling times are needed in the future to decide 
whether the shear modulus is continuous or discontinuous at 
the glass transition. 



and d — 2 dimensions, which we have sampled by means 
of MD and MC simulations, we have addressed the gen- 
eral question of whether the stress fluctuation method 
may be used through a solid-liquid transition where a 
well-defined displacement field ceases to be defined. De- 
liberately setting the experimentally motivated linear re- 
gression Eq. © as the fundamental definition, the shear 
modulus G has been computed from the strain and stress 
fluctuations (assuming a flxed finite measurement time 
window t) in ensembles where either a shear strain 7 = 
(NV7T- and NP7T-ensembles) or a (mean) shear stress 
T = (NVrT- and NPrT-ensembles) have been imposed 
(Fig. [T|). Working at constant mean normal pressure we 
have computed and compared the temperature depen- 
dence for various simple averages and fluctuations con- 
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tributing to the shear modulus in different ensembles. 
For the KA model only the currently available results for 
NV7T- and NP7T-ensembles have been presented. 

As has been stressed in Sec. Ill Bl and Sec. Ill CI the 
stress fluctuation representation Gtt at fixed shear strain 
does not rely in principle on a solid-like reference state 
for a displacement field of tagged particles if thermody- 
namics can be assumed to hold. It thus holds formally 
through the glass-transition temperature Tg up to the 
liquid regime (albeit with a trivial value G = 0) as does 
the better known relation for the compression modulus 
K (Appendix |X| 8]. As emphasized in Sec. Ill Dl impul- 
sive corrections due to the truncation of the pair poten- 
tials have to be properly taken into account, especially at 
high temperatures (Fig. U) , for the precise determination 
of the affine shear elasticity /xa = /^b ~ -Pox m- Con- 
firming Eq. p^ . it has been shown for the pLJ model 
that Grr at 7 = has the same long-time behavior as 
the conceptually more direct observables G^r and G^-y 
obtained using the strain fiuctuations at constant mean 
shear stress t — Q. The standard thermodynamic rela- 
tions comparing different ensembles appear thus to hold 
also in practice (at least within our numerical precision) 
through the glass transition for all T and this albeit 

1. some degrees of freedom get quenched at low tem- 
peratures on the time window t computationally 
accessible, 

2. all operational definitions of G are transient in the 
sense that they vanish for t — 00 (Figs. [T^ and [T5)) 
and 

3. it is not self-evident that in the latter limit and for 
large temperatures 7 and r can still be treated as 
a pair of conjugated thermostatistical variables. 

The latter point compactly epitomized by the thermody- 
namic relation Eq. ([8]) has been shown to hold, however, 
with remarkable precision up to very high temperatures 
for the pLJ model. 

As predicted by general Legendre transformation, G^-r 
vanishes for all temperatures T, if r rather than 7 is 
imposed (Figs. [SJ [71 and [TU)) . The shear-stress fiuctua- 
tions /^fIt- are thus given by the affine response under 
an external load /iA (Fig- HI, i-e. MfIt- reduces to a sim- 
ple two-point pair correlation function if pair potentials 
are considered as discussed in Sec. Ill CI The same holds 
above Tg for constant 7, since the boundary conditions 
are irrelevant for the liquid state. This symmetry with 
respect to the boundary conditions (r o 7) is broken at 
the glass transition for a large, but finite t: at constant 
7 the stress fiuctuations reveal a strong non-monotonous 
behavior with a clear maximum at Tg (Fig. 19]) . The shear 
modulus is the order parameter characterizing this sym- 
metry breaking. Alternatively, as we have discussed in 
Sec. Ill CI Gtt/ma = 1 — Mpl-y/MA may be seen as an 
order parameter comparing the ratio of the non-affine to 
the affine shear responses. Since /za > MfL for T < Tg, 



this implies that the stress fluctuations contain higher- 
order correlations reducing the free energy of the system. 

The increase of G below Tg is reasonably fitted for 
both models by a continuous cusp singularity, Eq. ([69l) . 
in qualitative agreement with recent replica calculations 
[H, . A jump discontinuity, as suggested by mode- 
coupling theory [Sll . l32j and another replica theory [s^ , 
is not compatible with the currently available data shown 
in Fig. [To] and Fig. ITT] However, as shown in Fig. [12] and 
Fig. [ini our data depend strongly on sampling time t and 
the computation of larger t could lead to a sharper tran- 
sition. Thus a jump discontinuity cannot be ruled out 
completely. (We are not aware of any other investigation 
of the sampling time effect for the shear modulus.) At 
present we believe, however, that it would be surprising 
if our data could be reconciled with a discontinuity at 
Tg. It seems more likely that below Tg, where the choice 
of the boundary conditions does matter as shown, the 
definition of the glassy shear modulus used in Refs. 
32] do not correspond to the key operational definition 
G^T- of the shear modulus considered by us. In any case, 
our work shows that any theory and numerical scheme 
used to determine the shear modulus using a stress or 
strain fluctuation relation (in the low-wavevector limit 
or at a small finite q) must cleary specify and take into 
account whether the shear stress or the shear strain are 
macroscopically imposed. Otherwise, such an approach 
is void. 



B. Outlook 

Beyond the presented work. Future work should 
clearly focus on the more detailed description of the sam- 
pling time dependence. By extrapolating appropriately 
for the large-i asymptotic behavior, this may allow to 
settle the theoretical debate. As emphasized above, one 
shortcoming of the present work is that the time-scale 
used in our MC simulations was slightly arbitrary, es- 
pecially if additional non-local particle moves are used 
in the low-T limit and if box volume and shape chang- 
ing strains are included. The comparison of i-dependent 
properties for different temperatures T becomes thus del- 
icate. Our MC study for the 2D soft particles should 
thus in any case be recomputed with Langevin thermo- 
stat MD dynamics as used for the data of the KA model 
presented. The t-dependence of the latter model, which 
is a better glass former than the 2D pL J model ^] , has 
still to be worked out for the key operational definition 
Gyr in the NVrT- or NPrT-ensembles |57!|. We plan also 
in the near future to reconsider more carefully our previ- 
ous investigation of the glass transition of polymer melts 
7] taking properly into account the impulsive truncation 
corrections and comparing the shear moduli G-yr ~ Gj-y 
(NPrT-ensemble) and Grr (NPrT-ensemble) around Tg. 

In addition, it should be rewarding to compare our 
shear moduli with the values obtained from the dis- 
placement of particles following the procedure chosen in 



19 



(a) 





<''•« y V' 








n. — 






\ \ \\ 


»y y;.i<:'.ji.''> 










•" T. / 


>- 



harmonic springs 
fluffy layer 

rigid rods moving only horizontally 
using transverse plane waves 

stiff layer 

beads glued on rods 

relaxation: horizontal plane waves 



(b) 
G(t) 



X 

T=const 



fixed bottom wall 

(c) t=const 
no or weak relaxation 

Oo o 

intermediate relaxation GC") " o 
^ frequency 



rapid stuess relaxation 



o o 00» 



log(t) 



FIG. 14; Outlook: (a) A simplified scalar model for per- 
manent and transient networks where beads (filled spheres) 
are glued transiently on rigid rods which are only allowed to 
move horizontally [5g |. (b) Shear modulus G{t) as a function 
of sampling time t at fixed temperature T for three different 
mean barrier heights {Eh)- (c) Possible shape of G{T) vs. 
temperature T for a fixed sampling time t. 



Ref. [3l| . This procedure relies on the idea that the par- 
ticles fluctuate around a well-defined reference position. 
In the low-T limit where the interaction network can be 
replaced by the dynamical matrix, this should yield the 
same results. However, for larger temperatures where 
the harmonic approximation must break down, this ap- 
proach is questionable (both for the analysis of experi- 
mental and computational data) and should be compared 
to the moduli obtained directly from the strain and stress 
fluctuations of the overall simulation box. 

A simple scalar model. We are currently working on 
an extremely simplified scalar model for permanent and 
transient networks in c? = 2 dimensions which may help 
to clarify the theoretical debate [58|. As sketched in the 
panel (a) of Fig. [TH in this model beads (filled spheres) 
are glued on rigid rods which are only allowed to move 
horizontally. In its most simple implementation this is 
done by means of global MC moves using transverse plane 
waves with a wavevector pointing in the y-direction. The 
glued beads are connected by ideal harmonic springs 
where the spring constants Ki and reference lengths Ri 
are chosen randomly according to fixed, narrow or broad 
distributions. The stiffness of the contacts between rods 
may differ (even including negative Ki and Ri as for the 
dynamical matrix studied in Sec. IIVB|) . The shear mod- 
ulus may be computed either from the shear strain fluc- 
tuations in the NVrT- or the shear stress fluctuations in 
the NV7T-ensemble. The shear stresses can be relaxed 
by either allowing the particles to move along the rods 
using longitudinal plane waves in the x-direction and/or 
by breaking and reconnection of springs. The latter step 
is currently done by applying a chemical potential for 



the springs. Physically, the study of such transient net- 
works is motivated by systems of hyperbranched polymer 
chains with sticky end-groups [35| or telechelic polymers 
in water-oil emulsions [3^, |33 . The associated relaxation 
frequency w; oc exp(— i?b/T) for a given spring I is as- 
sumed to be set by an activation barrier E\j which may 
differ for each spring. Depending on the distribution of 
this activation barrier the shear modulus G{t) thus de- 
cays more or less rapidly with sampling time t as sketched 
in panel (b) of Fig. [TH The shear modulus becomes con- 
stant, if the random percolating network is permanent 
(a;; = 0, Eb/T — )• oo) or the relaxation negligible (top 
solid line) . The system behaves as a liquid (bottom line) 
for large (mean) (w/), while a shoulder is seen in the 
intermediate frequency range. While this is conflrmed 
by preliminary simulations, we do at present not under- 
stand the role of the quenched noise injected in the model 
which is seen to strongly influence the shear stress fluc- 
tuations fiF (especially if fluffy layers happen to appear). 
Also the understanding of flnite system-size effects is of 
course crucial for such a low-dimensional model. (Only 
systems containing N = particles have been consid- 
ered so far.) As sketched in panel (c) the key question 
is to describe accurately for asymptotically long simula- 
tion runs and asymptotically large box sizes the shape of 
G{T) around the glass transition temperature Tg{{Eb)) 
where the shear modulus vanishes and /ip]^ should again 
have a maximum. 
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Appendix A: Compression modulus 

1. Operational definitions 

We have focused in the main part of this paper on 
the shear modulus G of our glass-forming model systems 
and the comparison of the three observables G-yr, G-yy 
and Grr- Being isotropic the complete macroscopic elas- 
tic description of our systems only requires the numeri- 
cal determination of one additional elastic modulus, the 
isothermal compression modulus 0, HO, 113 ; 



K 



-V 



dP 
dV 



dP 



(Al) 



characterizing the elastic response with respect to a 
volumetric (dilatational) strain fluctuation as shown in 
Fig.fTfd). As discussed in the standard textbooks [3.l40|. 
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the extensive variable is here the volume X — V ^ {V) 
and the conjugated intensive variable the negative pres- 
sure / = — P = — (-P) where V and P stand for the 
instantaneous values of the volume and the pressure. 
With these notations Eq. (jAip is obviously consistent 
with Eq. (IT2]) or Eq. (IT4|) . In analogy with G^r for the 
shear modulus G, the compression modulus K may be 
determined in an NP7T- or an NPrT-ensemble using the 
linear regression relation 



K^p = -V (^SVSp'j I (bV- 



(A2) 



where the index vj) indicates that both the measured in- 
stantaneous volume V and pressure P have been used. 
The corresponding correlation coefficient is 



(A3) 



Since from Eq. (Hi]) we have {WbP) = -k^^T if P is 
imposed, this implies that should be identical to the 
better known relation [401 



= k^TV/ (SV^-' 



(A4) 



which corresponds to the operational definition G^-^ dis- 
cussed above. The compression modulus K may of course 
also be determined using the pressure fluctuations in sys- 
tems of imposed volume V as sketched in Fig. [He). As 
(to our knowledge) first derived by Rowlinson [8j], the 
corresponding stress fluctuation formula reads 



Kpp = rjA ^ Vf with rjp = f3V ( SP 



(A5) 



being the fluctuation of the total normal pressure (corre- 
sponding to fip) and rjA the "affine dilatational elasticity" 
(corresponding to /xa). We shall give below a proper op- 
erational definition for 77A. Here we only note that due 
to the general Legendre transformation Eq. (|18p for in- 
tensive variables one knows immediately that 



K. 



This implies by comparison with Eq. (jA5[) that 



?7A = r/Flp 



(A6) 



(A7) 



in analogy to Eq. ((2T|) for the shear, i.e. rj^ can be di- 
rectly be determined from the pressure fluctuations in 
NP7T- and NPrT-ensembles. As before for the corre- 
lation coefBcient c-yi-, Eq. P^ . it follows from Eq. (jA7[) 
that 



(A8) 



Since c^p < 1, the affine dilatational elasticity sets an 
upper bound to the compression modulus K. 



Rowlinson's stress fluctuation formula 



Introduction. We reformulate here briefly the deriva- 
tion of the stress fluctuation relation Kpp for the compres- 
sion modulus K of an isotropic system at imposed volume 
V by Rowlinson [8|. With V{0) being the volume of the 
unperturbed simulation box we assume a small relative 
volume change 



V{e)/V{0) - 1 



(A9) 



with all coordinates equally deformed as, e.g. 

a;(0)^x(e) =x(0)(l-|-e)i/'^ (AlO) 

for the x-coordinate in d-dimensions. (The argument 
(0) denotes again the reference system.) Using a similar 
rescaling trick as for the shear modulus G in Sec. IIIBl 
the interaction energy Us{e) of a configuration s of the 
strained system is expressed in terms of the coordinates 
(state) of the unperturbed system and the explicit metric 
parameter e. Due to the volume change we have now also 
to take into account the ideal gas (kinetic) contributions. 

General excess contributions. We derive first the ex- 
cess contribution Pcx to the total pressure P = Pid -I- Pex 
and the excess contribution i^ex to the total compression 
modulus K = Kid + ^^^ex valid for an arbitrary conserva- 
tive potential. The general relations Eqs. (|26ll29p stated 
above for an imposed shear strain 7 still hold with the 
relative volume change e replacing 7. Using Eq. (|28p it 
follows for the excess pressure that 



p — _ 



9Pcx(U) 

dV 



1 aPex(e) 



U(0) 9(1 
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du'M 



U(0) de 



(All) 



defining the instantaneous excess pressure 42|. In the 
second step we have taken the limit e — )■ and have 
dropped the argument (0). The average is again eval- 
uated using the weights of the unperturbed system, 
Eq. dni). Using Eq. (gT]) and Eq. ^ one obtains for 
the compression modulus 



K,: 



V 



a2Pex(U) 1 52Pex(e) 



= {U'Jie))/V 



y(0)9(l + e)2 



13V (SP, 



(A12) 



which is again understood to be taken at e = 0. The 
first term 7yA,ex = {U'J{e)) /V in Eq. (jAl2p corresponds 
to the change of the system energy assuming an affine 
strain transformation for all particle positions. It gives 
the excess contribution 77a,cx to the total affine dilata- 
tional elasticity 77A = Va.m + ??A,ex- The second term 
corresponds to the excess contribution 77F.CX to the total 
normal pressure fluctuation rjF = ??F,id + ??f,gx mentioned 
in Sec. lA II It corrects the overprediction of the affine 
strain contribution 77a,cx |18|. 
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Ideal gas contribution. We remind that for an ideal 
gas the isothermal compression modulus A'id is given by 
the ideal gas pressure Pjd = k^Tp and that the ideal gas 
pressure fluctuation rjF,id = I3V{SP^^) becomes 77F,id|v — 
Pid at constant volume [5] . Interestingly, using again the 
general Legendre transform Eq. (|18l) one sees that 



?7F,id| 



?7F,id| 



id 



2Pid 



(A13) 



for the ideal gas pressure fluctuations at imposed pres- 
sure. 

Due the assumed additivity of the total system Hamil- 
tonian, the stress fluctuation formula for the total com- 
pression modulus may now be written 



Pid 



VA.. 



(A14) 



While this result is useful for constant- ensembles, es- 
pecially if MC simulations are used, it is yet not in a 
form exactly equivalent to the general Legendre trans- 
form Eq. (jA5|) . Since the ideal and excess pressure fluc- 
tuations are decoupled, {SPidSPcx) = 0, this implies 
Vf = VF,id + Vf,cx- One may thus rewrite Eq. (|A14[) 
as Kpp ~ rjx — rjF with 



2Pi 



id 



?/A,cx = ?7F,id|p + ?7A,ex 



(A15) 



where we have used Eq. (jAlSP in the second step. As- 
suming Eq. (jA15l) assures that rj^ becomes equivalent to 
?7fIp- This has the nice feature that not only the func- 
tional Kpp, Eq. (lASp . should give the correct compression 
modulus K in the NV7T and the NVrT ensembles, but 
also that Kpp must vanish properly if the "wrong" NP7T 
and NP7T ensembles are used, just as Grr vanishes for 
NVrT- or NPrT-ensembles. 

Pair interactions. Up to now we have not taken ad- 
vantage of the fact that the conservative interaction po- 
tential is assumed to be a pair potential, Eq. (|4T|) . This 
will be used now to express Pex and TyA.ex in terms of sim- 
ple (ensemble independent) two-point correlation func- 
tions. Equation (|A10I) implies that the squared distance 

between two particles transforms as 



r2(0)^r2(e)=r2(0)(l + e) 



2/d 



It follows that 



d2r2(e) 2(2 -d) 2 



de2 



d^ 



(A16) 

(A17) 
(A18) 



where we have taken finally for both derivatives the limit 
e — > 0. For the first two derivatives of a general function 
/(r(e)) with respect to e this implies 



oe d 



(9e2 



d2 



(A19) 



(A20) 



for e — > and dropping the argument (0) on the right- 
hand sides. According to the general result Eq. (|A11[) for 
the normal pressure and using Eq. (IA19P one confirms for 
pair interactions the expected virial relation [1, ^ 



(A21) 



for the instantaneous excess pressure where the interac- 
tions between pairs of particles i < j are again labeled 
by an index I. Using Eq. (jAlSP the excess contribution 
to the affine dilatational elasticity becomes 

?7A,cx = Vb + Pcx where (A22) 
^ (y^^ '"N''(n) + riu'jri) ^ (A23) 

is sometimes called the "hypervirial" contribution to the 
compression modulus fs', 'o'l . It corresponds to the first 
term indicated in Eq. (jA20l) . By comparing Eq. (IA23P 
with Eq. (|Tf| it is readily seen that for isotropic systems 

77B = (^1 + MB - ^Pox, (A24) 

i.e. T^B and /ib and, hence, the affine dilatational elastic- 
ity ?7a and affine shear elasticity ma are closely related. 
(This justifies to call 77B a "Born coefficient".) Please 
note that Eq. (R5|) or Eq. (|JT4l) together with Eq. (IX22|) 
are perfectly consistent with the relations stated in the 
literature p-H, HH H^. We emphasize finally that, as 
discussed for the affine shear elasticity ma in Sec. Ill Bl it 
is inconsistent to neglect the explicit excess pressure con- 
tribution to r]A in Eq. (jA22p but to keep the underlined 
term of the Born coefficient ?]b, Eq. (IA23P j60| . 



3. Numerical findings 



Affine dilatational elasticity. The affine dilatational 
elasticity t/a = '7a, id + ?7A,ex = 2Pid + Pox + Vb obtained 
for both models is presented in the inset of Fig. [13] as 
a function of the reduced temperature x — T/Tg. (The 
values may also be found in Table U and Table |lll) We 
remind that for determining the Born coefficient 77B im- 
pulsive corrections need to be considered as discussed at 
the end of Sec. IIIDI As shown by the small filled trian- 
gles, TjA compares well with the total pressure fluctua- 
tion ?7f|p = 2Pid -I- r?F,ex|p in an NPrT-ensemble. The 



observed data collapse is expected from Eq. (|A7p and 
Eq. (jAlSp . As indicated by bold and dashed lines for, 
respectively, the low- and high-T limit 77A decreases es- 
sentially linearly with temperature. The coefficients of 
the lines (not shown) are consistent with Eq. (|A24p re- 
lating 7] A and fiA and the fits given in Fig. [S] for fiA- We 
emphasize that the affine dilatational elasticity rjA (and 
thus Ma) is not only larger for the KA model but, more 
importantly that, is also increases more rapidly with de- 
creasing T than the pLJ model. 
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FIG. 15: AfSne dilatational elasticity riA{T) and compres- 
sion modulus K{T) for both models. Inset: rjA vs. reduced 
temperature x = T/T^. The small triangles indicate j^f for 
the pLJ model using the NPrT-ensemble confirming that 
Eq. (|A7|) holds. Main panel: The Rowlinson formula Kpp, 
Eq. (|Xn| . for the NV7T-ensemble is given for both the KA 
model (spheres) and the pLJ model (squares). For the pLJ 
model we indicate in addition K^p (crosses), K^^ (large trian- 
gles) and the rescaled correlation coefficient c^p (stars) which 
collapses on the other pLJ data as suggested by Eq. (|A8p. 



Compression modulus. The main panel of Fig. \W\ 
shows the best long-time values of the compression mod- 
ulus K for both models. The vertical axis is rescaled by 
77a. The rescaled stress fluctuation formula, Eq. (jA14p . 
for the NV7T-ensemble is given for both models. For 
the pLJ model we also indicate Kyp (crosses) and K^y 
(large triangles) for the NPrT-ensemble. No rescaling 
with 77a is necessary for the squared correlation coeffi- 
cient clp (stars) also presented. A perfect collapse of all 
pLJ data is found confirming the equivalence of all op- 
erational definitions used. For both models we observe 
a strong increase of K/rjx around a: ~ 1 with decreas- 
ing temperature. Interestingly, the rescaled data of both 
models is similar but not identical. The mentioned in- 
crease around a; sa 1 is clearly more pronounced for the 
KA model. Note that if the unsealed K is directly plotted 
against T, this effect becomes even stronger (not shown) 
consistently with the scaling of 77A presented in the inset. 
We note finally that for very low temperatures we observe 
for both models that the ratio Kj^f^ = c^^ approaches 
unity, i.e. the compression modulus is essentially dom- 
inated by the affine response. We remind that this is 
different for the reduced shear modulus for both models 
which is observed to approach G/ziA ~ 1/2 in the low-T 
limit (Fig. [TlT) . i.e. the shear stress fluctuations ^fL are 



more relevant than the pressure fluctuations 77F | y 
limit. 



this 



Sampling time dependence. For the pLJ model we 
present in Fig. [12] the sampling time dependence of var- 
ious operation definitions of K for the liquid limit at 
temperature T — 1.0. The horizontal axis indicates the 
sampling time t in units of MCS, the vertical axis is made 
dimensionless using ?7a as natural scale. As shown for 
Kyp and Kyy obtained in the NPrT-ensemble and for 
Kpp obtained using Eq. (UT4| in the NV7T-ensemble, 
all operational definitions applied in their natural ensem- 
ble converge rapidly within t ~ 10'^ MCS to the same 
plateau. (This is a general finding for all temperatures.) 
The plateau value for T — 1.0 is indicated by the bold 
solid line. Note that K/r]A is smaller unity as expected 
from the discussion at the end of Sec. lA 11 

The small filled symbols refer to Kpp computed in 
the "wrong" NPrT-ensemble. The small diamonds have 
been computed using directly Eq. (jA14l) . As indicated by 
the dashed line, this yields Kpp 2Pid for large times. 
Only if we compute Kpp = rjA — rjp using Eq. (jAlSI) . 
it does vanish properly as shown by the small triangles. 
This shows numerically that Eq. (|A14p . albeit perfectly 
fine for constant-!^ ensembles, is just not the correct for- 
mulation allowing to make manifest the general Legendre 
transformation Eq. (jA6p . This is, however, achieved us- 
ing Eq. (|A15p . This finding holds rigorously for all tem- 
peratures as can be seen from the scaling of the small 
triangles in the inset of Fig. [15] 
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FIG. 16: Rescaled compression modulus K(t)/r]A vs. sam- 
pling time t in MCS in the high-T liquid regime for pLJ par- 
ticles at P = 2 and T — 1 comparing K^p, Kw and Kpp in 
different ensembles as indicated. The filled symbols refer to 
Kpp obtained in the "wrong" NPrT-ensemble. 
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